Hopping in a Supercooled Lennard-Jones Liquid: Metabasins, Waiting Time 

Distribution, and Diffusion 
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We investigate the jump motion among potential energy minima of a Lennard-Jones model glass 
former by extensive computer simulation. From the time series of minima energies, it becomes clear 
that the energy landscape is organized in superstructures, called metabasins. We show that diffusion 
can be pictured as a random walk among metabasins, and that the whole temperature dependence 
resides in the distribution of waiting times. The waiting time distribution exhibits algebraic decays: 
r -1//2 for very short times and r _Q for longer times, where a « 2 near T c . We demonstrate that 
solely the waiting times in the very stable basins account for the temperature dependence of the 
diffusion constant. 



The energy landscape picture proposed more than 
thirty years ago by Goldstein has turned out to be a 
fruitful way of describing the complicated many-particle 
effects in disordered systems |2|. Starting from the joint 
potential energy landscape (PEL) V(x) of N particles 
as a function of their configuration x = {xi, ...xjy} one 
expects that the properties of the system at sufficiently 
low temperatures will be dominated by long residences 
near local minima of V(x) (inherent structures) with 
rare hopping events between them |3|. Recently it be- 
came clear for a model glass former that the strict hop- 
ping picture approximately holds for T < T c (landscape- 
dominated regime) |J,|5j,|y| where T c is the mode-coupling 
temperature |2j. However, even for higher temperatures 
T < 2T C many dynamic properties are still related to the 
properties of inherent structures (landscape-influenced 
regime) H 0, 

Thus, in both temperature intervals 
an observable like the diffusion constant D(T) should de- 
pend on the topography of inherent structures (IS) or, 
more generally, of basins (a basin of an inherent struc- 
ture is defined as the set of configurations that reach this 
minimum via steepest descent |llj|). Such an understand- 
ing is indispensable to grasp the underlying physics of the 
Adam-Gibbs relation [Hill. 

A simplified picture of glassy dynamics has been ex- 
pressed in phenomenological models in 0, 0, based 
on spatially uncorrelated hopping processes (random- 
walk) in configuration space. Then the whole temper- 
ature dependence is contained in the average waiting 
time (t(T)). The true dynamics of glass forming sys- 
tems, however, is expected to be more complicated. For 
example it is known that back- and forth correlations 
cannot be neglected and that the elementary jump dis- 
tances depend on temperature 0, ^^|- I n general, in 
a hopping approach the temperature dependence of the 
diffusion constant may be related to spatial and temporal 
aspects, as expressed by the relation 



D(T) 



a 2 (T) 
6N(t(T)Y 
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With this ansatz, we anticipate the important role of 
the mean waiting time and collect the spatial details of 
hopping in an effective jump width a(T). The latter 
involves (i) the average jump distance, (ii) correlations 
of jump widths with waiting times, and (iii) directional 
correlations of successive jumps. To our knowledge this 
decomposition into spatial and temporal contributions 
has not been systematically implemented within the PEL 
framework so far. A priori it is not clear to which de- 
gree the temperature dependence of a(T) is relevant; see, 
e -g-, [lil l2(j|. Some information about the waiting time 
distribution (WTD) has already been gained from the 
analysis of hopping processes of single particles in real 
space via computer simulations [2ll l22l|. In contrast, we 
consider hopping in configuration space, with the advan- 
tage of incorporating the full many-particle effects |2jj . 

In this paper, we present detailed information about 
the spatial and temporal aspects of hopping in a 
model glass former, and individually determine a(T) 
and (t(T)). We demonstrate that (i) only (r(T)) de- 
pends on temperature, (ii) hopping among single basins 
is not a random walk, whereas hopping among super- 
structures of minima (metabasins) is close to a random 
walk, (iii) (r(T)) is dominated by the long waiting times, 
due to a slow (approximately algebraic) decay of WTDs. 

In the present work, we investigate a binary mixture 
of Lennard-Jones particles (BMLJ), as recently treated 
by two groups ja, |24j ; see also |2j| . It is characterized 
by the interaction potentials V a p(r) — 4,e a /3[(a a p/r) 12 — 
(crap/r) 6 ] with the parameter set N = Na + Nb — 52 + 
13 = 65, a ab — 0.8ctaa, &bb = 0.88ctaa, £ab = 1-§£aa, 
(bb = 0-5tAA, r c = 1.8ajiA. Linear functions were added 
to the potentials to ensure continuous forces and ener- 
gies at the cutoff r c . Units of length, mass, energy, and 
time are ctaa, m , £aa, and y/ma AA /eAA, respectively. 
For convenience, though, we will omit units here. We 
use Langevin molecular dynamics simulations (MD) with 
fixed step size, A 2 = 0.015 2 = 2kQTAt/m(, equal parti- 
cle masses m, friction constant £, and periodic boundary 
conditions at a density of p = 1.2. The friction constant 
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C = 2/0.015 2 is chosen so that At = l/T. Due to the dif- 
ferent type of dynamics, the absolute values of times and 
diffusion constants are different from those found within 
Newtonian dynamics simulations. The mode-coupling 
temperature is T c = 0.45 ± 0.01 in this model system 
(compare |25jl. 

For the analysis of dynamics in configuration space it 
is essential to use small systems because otherwise many 
interesting effects will be averaged out [lflj. The rele- 
vance of small systems has been also pointed out by other 
groups; see, e.g. [2^, [2^. On the other hand, the system 
should not be too small in order to avoid major finite 
size effects. N w 60 turns out to be a very good compro- 
mise for binary Lennard- Jones mixtures, whereas N < 40 
already displays major finite-size effects j2?|. Here we 
choose TV = 65. To back those findings,we have carried 
out an extensive study of finite-size effects for systems 
of N = 65, 130, and N = 1000 particles 28]. It turned 
out that the N = 65 system is nearly identical to the 
bulk (N = 1000) above T c . Since well-equilibrated runs 
of N > 130 are lacking below T c , finite-size effects can- 
not be excluded there at the present stage. However, 
this does not affect the main results of this paper. More 
important is the question of good equilibration at each 
temperature. Have the runs been long enough to sam- 
ple the PEL sufficiently? Above T c this is uncritical, 
which can be seen from the fact that each run comprised 
at least 850 a-relaxation times. A more detailed check, 
involving the lifetimes and distribution of metabasins, in- 
dicates that runs down to T = 0.435 are feasible with the 
available computer power. 

By regular quenching the MD trajectory x(t) to the 
bottom of the basins visited at time t, as proposed by 
Stillinger and Weber, we obtain a discontinuous tra- 
jectory £(t). In this way, one discards the more or 
less complicated vibrational part x(t) — £(t) of motion, 
only keeping the visited minima as 'milestones'. The 
one-particle diffusion constant can be also determined 
from the squared displacement of inherent structures via 
D = lim t ^ 00 <(£(*) - m?) /GNt. 

How to resolve the elementary hopping events? Since 
computer time prohibits to calculate for every time 
step, we normally find ourselves in the situation of hav- 
ing equidistant quenched configurations with, say, 
U+i — i< ~ 10 5 MD steps. If the same minimum is found 
for times U and tj , we need not care about transitions in 
the meantime, because no relaxation has occurred there. 
If, in contrast, f (t$) ^ we must not expect £(tt+i) 
to be the direct successor of £(U), since many other min- 
ima could have been visited between ti and ti+j.. There- 
fore, further minimizations in this time interval are neces- 
sary. For reasons of efficiency, we apply a straightforward 
interval bisection method, which locates all relevant tran- 
sitions with an accuracy of 1 MD step. Although compu- 
tationally demanding, this has proven most efficient for 
resolving the relevant details of hopping on the PEL. 
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FIG. 1: The time series of minima energies measured for 
the BMLJ system of TV = 65 particles, T = 0.435. The 
distance between minimizations before interval bisection is 
10 5 MD Steps. The length of the total run is 2 x 10 9 MD Steps. 
Top: time window covering a quarter of the total run. Bot- 
tom: magnification by a factor of fifty. 



As demonstrated in |10j, the time series of potential 
energies e(t) = V (£(*)) reflects well the character of dy- 
namics in the supercooled state. For T = 0.435, e(t) is 
shown in Fig. ^ from which we note a remarkable struc- 
ture in e(t). The system is trapped in some stable config- 
urations for long times, during which a small number of 
minima is visited over and over again. Obviously, these 
minima form superstructures, which, following we de- 
note metabasins (MBs). One may imagine that minima 
of long-lived MBs are organized in funnel- like structures 
so that the system is stuck there for a long time. It has 
been argued that the occurrence of /3-relaxation at low 
temperatures is due to such substructure of the PEL |3] . 
This is supported by the real-space signature of MBs as 
reported by Middleton and Wales [2t| . Formally there is 
no unique way to define MBs for a given PEL due to the 
lack of a strict time scale separation. Here we take the 
pragmatic view and let the system decide by its MD run. 
The intuitive notion of MBs from Fig.Q]can be cast into 
an algorithm in a straightforward way, see [Toj |36l ] . One 
major advantage of analyzing MBs rather than basins is 
that the charming simplistic picture of a random walk 
in configuration space will be better fulfilled since direct 
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FIG. 2: Squared displacement (_R 2 (n)) after n MB jumps 
for the temperatures T = 0.435,0.466,0.5,0.6,0.8, and 2.0. 
Also shown are three curves (T = 0.5, 0.6, and 0.8) for jumps 
among single basins (lower curves). The latter data have been 
obtained by extremely frequent minimizations, which resolve 
nearly all IS transitions. We have included lines of slope 1. 



back- and forth correlations are already taken into ac- 
count. It turns out (see Fig. ^) that long-lived stable 
MBs are separated by bursts of rapid transitions among 
higher minima which look like fountains. They cannot be 
detected without interval bisection. The waiting times of 
the MBs range from a few MD steps to many millions of 
them. Moreover, comparing Fig. ^ ( a ) an d (b), we find 
a certain self-similarity in e(t) when inspected on differ- 
ent time scales. The distribution of MB lifetimes r will 
be denoted ip(r,T), its first moment (r(T)) being a key 
quantity for the following considerations. 

Our goal is to find an expression for the effective jump 
width a(T) of Eq. QJ. To this end, as a generalization 
of we define the MB inherent structure £ MB (i) as 
the mean of £(/;) over the MB lifetime. The key idea 
now is to introduce the squared distance (i? 2 (n)) after 
n jumps |30| . averaged over different realizations. This 
quantity is purely spatial since it does not involve any 
time scale. In the limit of large n one obtains due to the 
central limit theorem lim n ^ 00 (_R 2 (n))/(^ B (r7, (r))) = 1 
where (^(n (t))) is the average squared displacement 
after time n (r). Thus, 



D{T) = lim 



6Nn (t) 



lim 



6 77 



1 



iV(r) 



where the first factor may be identified as a 2 (T)/6. Note 
that both factors are independent of system size because 
IS transitions are localized ((i? 2 (l)) = 0(1)) and mean 
waiting times decrease with system size ((r) = 0(1/N)). 
We may now calculate a(T) from the simulations, re- 
sults are shown in Fig. The most important ob- 
servation is that a(T) is temperature independent for 
T < 1. Interestingly, this is not affected by the variation 



FIG. 3: Arrhenius plot of the inverse one-particle diffusion 
constant 1/D(T) and the mean waiting time (r(T)} multiplied 
by a constant (a 2 = 1.0). Error bars are of the order of the 
symbol size. 



of the elementary jump width (i? 2 (l)), which increases 
with temperature, a fact that was recently observed by 
Schulz et al. [l^. A possible explanation for this might 
be found in the increasing population of higher-order sta- 
tionary points []|, which are known to be more distant 
to neighboring minima, but also provoke many 'book- 
keeping' IS transitions, thus resulting in larger backward 
correlations [23| . 

In any event, the constancy of a(T) in the landscape- 
influenced regime T < 1 implies that the tempera- 
ture dependence of D(T) follows alone from (r(T)), i.e. 
D(T) oc 1/ (t(T)). This simple picture breaks down for 
T > 1 where the explored regions of configuration space 
probably have a completely different structure. It has 
to be noted that the constancy of a(T) for low T relies 
heavily on our resolution of all elementary IS transitions 
leading to relaxation. 

A further insight from Fig. [21 is that the dynamics on 
the level of MBs is basically a random walk except for mi- 
nor back-correlations for n < 5. As seen from the figure, 
these correlations are present between single basins, the 
consequence being a significant deviation from the rela- 
tion (R 2 (n)) oc n. Also note the oscillations in the single- 
basin (R 2 (n)^ at small n, which result from the back- 
and-forth motion within MBs. More importantly, the 
single-basin curves do not have the same large-n limit so 
that the effective jump length on the level of single basins 
would be temperature dependent. It remains unclear why 
for a very small LJ system (N = 32) correlations among 
adjacent basins are irrelevant for T m T c 23], and why 
intra- rather than inter-basin dynamics is deemed to be 
the key to the understanding of diffusion 31]. 

We can check the relation D oc 1/ (t) within our simu- 
lations. Figure 3 shows that it is indeed well fulfilled for 
T < 1 while for T = 2 we find the expected deviations. 
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PEL has provided new insights into the mechanism of 
diffusion in supercooled liquids. As we have seen, the 
emergence of long-lived MBs is the reason for the slow- 
ing down of molecular motion in our supercooled model 
liquid. As a next step the waiting time should be related 
to the respective MB energies to establish a connection 
between energy (thermodynamics) and dynamics in the 
spirit of the Adam-Gibbs relation (see the subsequent 
paper |35j). 

We gratefully acknowledge helpful discussions with 
J. P. Bouchaud, T. Odagaki, D.R. Reichman, R. Schilling, 
H.R. Schober, and H.W. Spiess. Funding was granted by 
the Sonderforschungsbereich 262. 



FIG. 4: Distributions of waiting times, tp(r,T), for T — 
1.0,0.8,0.6,0.5, and 0.435 (from left to right). Curves have 
been shifted to overlap for small r. Lines corresponding to 
algebraic decays with exponents 0.5 and 2 are shown as guides 
to the eye. The arrows mark t*(T), i.e. the waiting time 
with the property f^?, T s Artp{r)r = 0.9 (r). Inset: Power-law 
exponent a(T) from fits to the long-time decay. Within the 
possible accuracy, a(T) » T/T c + 1. 



In the remaining part of the paper, we discuss the prop- 
erties of the WTDs ip(r,T), see Figure 4. For short r, 
all curves exhibit a power-law behavior with exponent 

At r k 10 3 — 10 4 a crossover 



— 1/2, similarly to |J 
to the faster decay <p(r,T) oc t _q ( t ) can be observed. 
For T w 0.45, one finds a « 2.0 for which the expec- 
tation value (r) would diverge. However, the behavior 
ip(r) oc r~ a cannot extend to infinity. Due to the finite 
number of MBs in the system, there exists a maximum 
effective barrier i? max , giving rise to an exponential cutoff 
at some minimum rate 7 min . 

The slow decay of the WTDs leads to the important 
observation that the mean waiting time (r) is dominated 
by the contributions from large r (see arrows in figure 4). 
For example at T = 0.435, ninety percent of (r) are made 
up by the six percent longest waiting times, which are 
basins of lifetimes greater than 0.5 x 10 6 MD Steps. This 
may come as a surprise because (i) the short-lived MBs 
are much more numerous than the long-lived and (ii) 
one might intuitively think that the diffusion constant 
and thus (r) is dominated by the fast particles. The 
result is in qualitative agreement with the approach of 
Wolynes and Xia who regard the relaxation of long-lived 
local structures as the time-determining step 33]. 

Interestingly, the algebraic decay of the WTDs fol- 
lows for some theoretical models of diffusion with built-in 
traps. This is the case for Bouchaud's trap mod el l lol and 
the trapping diffusion model of Odagaki et al. 15]. A re- 
cent comparison of WTDs in the Lennard-Jones system 
with that of trap models can be found in [34| . 

In conclusion, the detailed analysis of hopping on the 
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